eV=[0.64 0.77 0.89 1.02 1.14 1.26 1.39 1.51 1.64 1.76 1.88 2.01 2.13 2.26 2.38 2.50 2.63 2.75 2.88 3.00 3.12 3.25 3.37 3.50 3.62 3.74 3.87 3.99 4.12 4.24 4.36 4.49 4.61 4.74 4.86 4.98 5.11 5.23 5.36 5.48 5.60 5.73 5.85 5.98 6.10 6.22 6.35 6.47 6.60]';
%gold
%n=[0.92 0.56 0.43 0.35 0.27 0.22 0.17 0.16 0.14 0.13 0.14 0.21 0.29 0.43 0.62 1.04 1.31 1.38 1.45 1.46 1.47 1.46 1.48 1.50 1.48 1.48 1.54 1.53 1.53 1.49 1.47 1.43 1.38 1.35 1.33 1.33 1.32 1.32 1.30 1.31 1.30 1.30 1.30 1.30 1.33 1.33 1.34 1.32 1.28]';
%k=[13.78 11.21 9.519 8.145 7.150 6.350 5.663 5.083 4.542 4.103 3.697 3.272 2.863 2.455 2.081 1.833 1.849 1.914 1.948 1.958 1.952 1.933 1.895 1.866 1.871 1.883 1.898 1.893 1.889 1.878 1.869 1.847 1.803 1.749 1.688 1.631 1.577 1.536 1.497 1.460 1.427 1.387 1.350 1.304 1.277 1.251 1.226 1.203 1.188]';
%silver
n=[0.24 0.15 0.13 0.09 0.04 0.04 0.04 0.04 0.03 0.04 0.05 0.06 0.05 0.06 0.05 0.05 0.05 0.04 0.04 0.05 0.05 0.05 0.07 0.10 0.14 0.17 0.81 1.13 1.34 1.89 1.41 1.41 1.38 1.35 1.38 1.31 1.30 1.28 1.28 1.26 1.25 1.22 1.20 1.18 1.15 1.14 1.12 1.10 1.07]';
k=[14.08 11.85 10.10 8.828 7.795 6.992 6.312 5.727 5.242 4.838 4.483 4.152 3.858 3.586 3.324 3.093 2.869 2.657 2.462 2.275 2.070 1.864 1.657 1.419 1.142 0.829 0.392 0.616 0.964 1.161 1.264 1.331 1.372 1.387 1.393 1.389 1.378 1.367 1.357 1.344 1.342 1.336 1.325 1.312 1.296 1.277 1.255 1.232 1.212]';
hbar=jcb_hbar;
%h=hbar*2*pi;
c=jcb_c;

q=1.602176487e-19;


%J=eV*q;
w=eV*(q/hbar);
lambda=2*pi*c./w;
lambda_nm=1e9*lambda;

w_interp=linspace(min(w),max(w),1000)';
lambda_interp=2*pi*c./w_interp;
lambda_nm_interp=1e9*lambda_interp;
n_interp=spline(w,n,w_interp);
k_interp=spline(w,k,w_interp);

% figure(1);
% plot(lambda_nm_interp,n_interp,'r-');
% hold on
% plot(lambda_nm,n,'bx');
% hold off
% title('n');

% figure(2);
% plot(lambda_nm_interp,k_interp,'r-');
% hold on
% plot(lambda_nm,k,'bx');
% hold off
% title('k');

epsilon=(n_interp+1i*k_interp).^2;
e2 = 1; %air
%e2 = 1.54385; % fused silica (quartz)
%e2 = 1.45810; % amorphous silica (glass);
wp=17.5e15; %silver
epsinf=8.5; %silver
epsilon_drude=epsinf-(wp^2)./(w_interp.^2);
%epsilon=epsilon_drude;

kx=sqrt(epsilon*e2./(epsilon+e2)).*(w_interp/c);
kx_drude=sqrt(epsilon_drude*e2./(epsilon_drude+e2)).*(w_interp/c);

[tmp maxind]=max(real(kx_drude));
[tmp minind]=min(real(kx_drude));



xmult=1e7;
ymult=1e15;
%figure(1)
clf;

hold on
%plot(imag(kx)/xmult,w_interp/ymult,'r-');
L(1)=plot(sqrt(e2)*([0 50]*ymult)/xmult/c,[0 50],'k-');
L(2)=plot(sqrt(1.45810)*([0 50]*ymult)/xmult/c,[0 50],'k--');
L(3)=plot(real(kx_drude(1:maxind))/xmult,w_interp(1:maxind)/ymult,'r');
plot(real(kx_drude(minind:end))/xmult,w_interp(minind:end)/ymult,'r');
% we set a linewidth which is !=1 so that it gets a width in the pdf
% output. However, we don't want it to be too thin, and  1 is actually a
% good width. Hahaha.
L(4)=plot([0; real(kx)/xmult],[0; w_interp/ymult],'b-','linewidth',1+eps);
ax(1)=gca;
leg=legend(L(1:2),{'light line in air','light line in glass'});
set(leg,'location','SouthEast','box','off')
ylim([0 10])
xlim([0 ceil(max(real(kx)/xmult))])
xlabel(['k_x \times 10^{'    num2str(log10(xmult)) '} / m^{-1}'])
ylabel(['\omega \times 10^{' num2str(log10(xmult)) '} / rads^{-1}']);

ax(2)=axes('color','none','box','off','xtick',[]);
ax(3)=axes('color','none','box','on','xtick',[],'ytick',[]);
set(ax,'position',get(ax(1),'position'));
set(ax(1:2),'box','off');
set(ax(2),'Yaxislocation','right',...
    'Ytick',w2lam(fliplr([200,250,300,400,600,1000,2000]*1e-9))/ymult,'ylim',get(ax(1),'ylim'));
set(ax(2),'YtickLabel',num2str(w2lam((get(ax(2),'Ytick')*ymult))'*1e9));
ylabel(ax(2),'\lambda_0 / nm');
set([ax(2) ax(1) get(ax(1),'xlabel') get(ax(2),'ylabel')...
    get(ax(1),'ylabel')],'FontName','Times New Roman','FontSize',10);
set(ax,'position',[0.07 0.145 0.84 0.83])


%%

% set([gcf ax(1) ax(2) ax(3)],'units','pixels');
%set([ax(1) ax(2) ax(3)],'position',[padding{1} [4 3]*scaling]);
figurename=fullfile(fileparts(mfilename('fullpath')),'DispersionRelation');

set(gcf,'color','white');
set(gcf,'PaperPositionMode','auto','units','centimeters',...
       'PaperSize',[16 10],'PaperPosition',[0 0 16 10]);
print('-dpdf','-r1000',figurename);
%print('-dpng','-r1000',figurename);
